General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



I 

I \ 


- 'it 




38>0 














r , 






.. s-r?L: 


. 




rv 2.^ ' > » 


e&ti 




a V 


% 


i 

1 


m 


r 












. 




and Space Administration 




X ; ,? ■ 






H O U S "LjCU-bL . TEXAS 


i K ■ 


= \ \ \ \/\ l I . 

u (SaSACR" or TMXOR ad NUMBER) 


(CAItOO* i I 


£<?ui?£ 


. 

■ 


Manned Spacecraft Center 


; Y’"/ : 

I 




" C .-'v .. 


Y - .. > 


v. 


I 




• -s- 




. 


* 

... 




* 

* 


V }-: . 





MSC-IN-64-ED-7 


MSC INTERNAL TECHNICAL NOTE 


A MODIFICATION OF HERRICK'S SOLUTION 
OF THE TWO-BODY PROBLEM FOR ALL CASES 


BY 


William H. Goodyear 
Lawrence F. Guseman, Jr. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
MANNED SPACECRAFT CENTER 
HOUSTON, TEXAS 
July 28, 1964 


MSC-IN-64-ED 


MSC INTERNAL TECHNICAL NOTE 


A MODIFICATION OF HERRICK’S SOLUTION 
OF THE TWO-BODY PROBLEM FOR ALL CASES 


Prepared by 

RTCC Project 

Prepared 

Theory and Analysis Office 




NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
MANNED SPACECRAFT CENTER 
HOUSTON, TEXAS 
July 28, 1964 


TABLE OP CONTENTS 


Page 

SUMMARY X 

LIST OF SYMBOLS 2 

INTRODUCTION 3 

SECTION 1: DERIVATION OF GENERAL SOLUTION PROM 

ELLIPTIC CASE 


1.1. Definition of the Transcendental 

Functions 6 

1.2. Properties of the Transcendental 

Functions 7 

1.3. Kepler’s Equation and Its Deriva- 
tives 9 

1.4. Closed Form Expressions for the 

Series Solution 12 

SECTION 2: A METHOD FOR COMPUTATION OF COORDINATES 

2.1. Initial Computations and Start 

of Iterations 15 

2.2. Evaluation of the Transcendental 

Functions 17 

2.3. The Solution of Kepler’s Equation 18 

2.4. Computation of Coordinates 21 

SECTION 3: A METHOD FOR COMPUTING PARTIAL DERIVATIVES 

3.1. Outline of Derivation of Partial 

Derivatives 23 

3.2. Evaluation of Parameters and 

Periodicity Computations 24 

3.3. Evaluation of the Four by Four 

Matrix 25 

3.4. Computation of the Partial 

Derivatives 26 

CONCLUDING REMARKS 28 


REFERENCES 


29 


SUMMARY 


A solution of the two-body problem given by Herrick 
for many types of orbits is modified by using a different 
independent variable. The modification changes Herrick’s 
transcendental functions and his form of Kepler’s equation. 
Th* end result is a more general solution of the two-body 
problem in that it anplies to both attractive and repulsive 
forces of any magnitude. 

An outline is given for deriving the general modified 
solution from the equations for an elliptic orbit. This re- 
quires definition of transcendental functions which are then 
used to express Kepler’s equation and give closed-form ex- 
pressions for the series solution to the differential equa- 
tions. A method is described in detail which outlines the 
computation necessary to determine coordinates at a given 
time from their known values at a given reference time. 
Formulas are also given for computing the partial deriva- 
tives of each of the resulting coordinates with respect to 
each of the reference coordinates. 


LIST OP SYMBOLS 


t 

x, y, z ^ x, y, z 


h 


w 

a 

E 

e 

r 

• 

r 

v 


S 0 , ^2* 

si, sj 


P 

t* g 
® • 
g 

At 


- time 

- rectangular position and velocity co- 
ordinates. When subscripted they 
refer to a particular point in time. 

- sum of the potential and kinetic 
energy 

- gravitational constant 

- semi -major axis 

- eccentric anomaly 

- eccentricity 

- magnitude of position vector 

- magnitude of change in position vector 

- magnitude of velocity vector 

- defined parameter used to generalize 
the equations for elliptic motion 

- defined transcendental functions 

- period of the orbit 

- power series in ( t — 1 0 ) 

- time derivatives of f and g . 

- residual used in solution of 
Kepler’s equation 
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INTRODUCTION 

A .number of proven methods exist for obtaining a space 
trajectory from given initial conditions at some reference 
time (Encke'St Cowell’s, Variation of Parameters, etc.). Most 
of these methods start with the undisturbed two-body space 
trajectory and consequently the need for a logically-simple , 
numericaliy-accurate solution of the two-body problem is 
apparent . 

The usual approach to solving the two-body problem is 
to develop a separate method of solution for handling each 
of the different types of orbits encountered. The choice of 
the method to be used is made by logically testing numerical 
values of certain orbit-defining parameters which have been 
computed from a set of Initial conditions. 

In practice it is sometimes difficult to determine numeri- 
cally the exact values of parameters for which one method is 
chosen as opposed to another (e.g., the problem of choosing 
the "best" when the energy is very small in absolute value, 
or zero). Consequently, it is advantageous to have a single 
method for solving the two-body problem which is continuous 
for all values of the orbit-defining parameters, thereby 
eliminating the logic associated with making a choice of 
methods. 

A solution of the two-body problem is given by Herrick 11 [1] 
for all cases in which the constant y in the differential 
equations 

* When this work was completed, the authors were not aware 
of the work done by K. Stumpff [2]. 


( 1 ) 


- l\ - 

x * -yx/r 3 

y * -yy/r 3 r = V x z + y* t z 2 

z * -yz/r 3 

Is positive and relatively large. His solution may be modi- 
fied to include all values of y by utilizing the parameter 

.(E-E ) 

* - .Q 

✓57a (2) 

rather than /a(E-E 0 ) to generalize the equations for ellio^- 
tic motion. 

For an elliptic orbit, E is the eccentric anomaly for 
any time t at which the position coordinates are x, y, z, 
and the velocity coordinates are x, y, z, and E 0 is the 
value of E for a particular reference time t 0 at which 
the position coordinates are x Q , y 0 , z 0 and the velocity 
coordinates are x 0 , y Q , z Q . Also, a is the semi-major 
axis of the orbital ellipse, and 

y/a » -2h (3) 

where the negative constant h is the sum of the kinetic 
and potential energy. 

h * v 2 /2 - y/r * Vg/2 - y/r Q ( 4 ) 

Here the square of the magnitude of the velocity is 

v 2 * x 2 + y 2 + z 2 (5) 

and r 0 , v 0 are values of r, v at the time t 0 . 


5 - 

The solution of the differential equations (1) ex- 
presses the relation between the coordinates x, y, e, x, 
y, z at arty time t and the coordinates x 0> y 0 , z 0 , x 0 , 
y 0 , z 0 at a reference time t Q . In Section 1, the equations 
of this solution for elliptic orbits are generalized to in- 
clude all types of orbits by using y rather than E or 
(E~E 0 ) as the variable in Kepler* s equation. This requires 
the expression of Herrick *s transcendental functions in 
terms of \p . In Section 2, a method is described for com- 
puting the coordinates x l9 y x , z^ x^, y x , i x at a given 

• • 

time t x from the known coordinates x 0> ( y 0 , x Q> y 0 , 
z 0 at a given reference time t 0 . In Section 3, formulas 
are given for computing the thirty-six partial derivatives 
of each of the coordinates Xj, y l# z , x^, y , z with 
respect to each of the coordinates x 0 , y 0 * z 0> x 0 , y 0 , z 0 . 
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1. DERIVATION OF GENERAL SOLUTION FROM ELLIPTIC CASE 
1.1. Definition of the Transcendental Functions 

The six transcendental functions S , S , S , S , 

0 12 3 

S , S are defined below In terms of narameters for 
** 5 

elliptic motion, but are expressed In terms of 2h 
and ^ for the general solution 


S h eos(E-E 0 ) 


(6a) 


- 1 - (E-E 0 ) 2 /2! + (E-E 0 )V'I! - (E-E 0 ) 6 /6! + ... 
<* 1 + (-^ 2 )/2! + <-HjP 2 > 2 /J< 1 + (-^ 2 ) 3 /6! + ... 

- 1 + (2h) 1 i|» 2 /2! + ( 2h) 2 <>V 1 l I + (2h)H 6 /6! + ... 


S = sln(E-En) » (E-E n ) sin (E-En ) 
% 


/y/a 


/w7a (E-E. ) 


(6b) 


- *[1 - (E-E 0 ) 2 /3! + (E-E 0 )V5! - (E-E 0 ) 6 /7I + ...] 

- ^C1 + (-£* 2 )/3! + (-^ 2 ) 2 /5! + (-^ 2 ) 3 /7t + ...] 

- ill 1 /!! + (2h) l Di 3 /3! + {2h) 2 K> 5 /5! + (2h) 3 + 7 /7 ! + ... 


S = 1 — COS (E-E n ) 

’ 2 ' ( /y7a ) 2 


S - 1 

-fl. 

^2h 


( 6c) 


i|« 2 /2 ! + (2h) 1 ! + (2h) 2 i|i 6 /6! + (2h) 3 i|i 9 /8! + ... 


S 


3 = 


(E-E fl ) - sln(E-E a ) _ Si-i|i 

Th “ 


(6d) 


(/y7a) 3 

i|i 3 /3 ! + (2h) 1 i|i s /5! + ( 2h) 2 i|i 7 /7 ! + (2h) 3 <J/ 9 /9! + ... 


_ (E-E ) 2 /2 1 - [ 1-cos (E-E )] S- -i|i 2 /2 ! 


(6e) 


(Ai/a) 


2h 


+ (2h) 1 i|i 6 /6 ! + (2h) 2 * 8 /a ! + (2h) 3 t ie /10! + 


• t • 



- 7 - 


g . (E-En)3/31 - C(E-E„) - Sln(E-E n ) ] _ S 3 -» 3 /3 1 (fif) 

5 ' </iT7iT) s 2 h 

- <!> S /5I + (?.h) l <h 7 nt + (2h) z i|iV9! + ( 2 h) 3 * ll /ll! + .. 
These transcendental functipns of ^ and h for 
the general solution replace the trigonometric func- 
tions of E or (E-E 0 ) that apply only to elliptic 
motion. 


1.2. Properties of the Transcendental Functions 

The transcendental functions Sj, S 2 , S 3 of h and 
ip have the following important properties: 


3 S, - s 

3 i)j 


aS i - ys 2 - ] .S 3 
3 h 

8S2 - S j 

3 ip 

(7a) 

* S 2 m *PS 3 - 2. Si* 
3h 

8S „« s, 

3\p 


9 ” 3 *S 5 


The use of these relations is the reason for defining 
the functions Sqj Sj and Si*, S 5 since the general 
solution of the differential equations ( 1 ) may be ex- 
pressed in terms of S 2 and S 3 alone. 

In the case of an elliptic orbit (for which h is 
negative) let 


where 


(t-t 0 ) * (t-to) - mP 


;7T 


2itp/(/ -2h) 3 


( 6 ) 




( 9 ) 


— 8 — 


is the period of the orbit, and m is a negative, 
zero, or positive integer which is chosen to mini- 
mize |(t-t 0 )| . Then 

( E -E 0 ) * (E-B 0 ) - 2irm (10) 

rather than (E-Eo) will he related to (tT-to), 
rat'her than (t-t 0 ), by the elliptic form (13) of 
Kepler’s equation* However, the coordinates x, y, 
z, x, t y, z are the same for (E-Eo) and (E-Eo) be- 
cause of the neriodicity of the orbit. 

Similarly, in the special case of an elliptic orbit, 

F ,B E ~ E n = * - m — (11) 

/iI7a / -2h 

rather than t will be related to (t-to), rather 

than (t-t 0 ), by the general form (16) of Kepler’s 

* * * 

equation. Also, the coordinates x, y, z, x, y, z 
are the same for ij# and iff because of the periodicity 
of art elliptic orbit. The definitions (6a) - (6f), 
(10), and (11) can then be used to show that 

50 * So 

51 = Si 

Sg * ^2 

5 3 * S 3 + m[2ir/(/^2h) 3 ] (12) 

5 4 » S 4 + m[ 2ir/(/-2h) 3 ] (i|/+l|/ )/2 

5 5 = S 5 + m[2Tr/(/-5h) 3 ][^ 2 +^+^2 )/6 + l/2h] , 



— 9 — 

where S Q , Sj, $ 2 , S 3 , S 4 , S g are respectively 
the transcendental functions (6a) through (6f) of 
5T and h rather than $ and h . 


1.3. Kepler ? s Equation and Its Derivatives 

Kepler’s equation for an elliptic orbit may be ex- 
pressed in the form 


(t-t 0 ) 


. M-M 0 

' w a 3 

= (E-e slnE) - (E 0 - eslnEp ) 
/ p/a 3 


(13) 


= a(l-eoosEo) + /5ae nln En l-cos(E-E 0 ) 

• / w7a ( /57a) 2 


+ pecosEo ^- 5 ^- 8 - 1 - ^ . ^ 

(/57a) 3 

In these equations, e is the eccentricity of the 
orbital ellipse. One expression for (M-Mq )// y/a 3 
is in terms of E and the other expression is in 
terms of (E-E 0 ). These two expressions are trigo- 
nometric identities, but the second is more general 
in that it is valid for circular and near-circular 
orbits . 

The parameters r, (rr) and (rv 2 - m) for an ellip- 
tic orbit may be expressed in the forms 

r = + /x z +y z +z z (l^a) 


a(l-ecosE ) 
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a(l - ecosEo) + /ya eslnEo ??- n 

/p7a 

1-cos (E-E fl ) 


+ pe cosE 0 


(/p7a) 2 


(14b) 

(rf) = xx + yy + zz 
* /piT eslnE 

» /pa e sinE 0 cos (E-E 0 ) + pe cosE ff 0 ^ 

/p7a 

(rv 2 - p) = + /x 2 + + z z (x 2 +y 2 +z 2 ) - p (14c) 

* pe cosEo 

« pe cosE cos(E-E ) + (-l±)/pa eslnE s ^ n ^ , T. J — 

0 0 a /57a 

Equation (14a) Is r times the time derivative of 

Kepler ? s equation (13). Similarly, (14b) is r times 

the time derivative of (l4a), and (14c) is r times 

the time derivative of (14b). For the particular 

case where t is to and E is thus Eo, equations 

(l4a) - (l4c) become 


r 0 = a( 1 

o 

PO 

m 

o 

0 

<u 

1 

(15a) 

(r o i 'o ) = 

/pa e sin E 0 

(15b) 

(r 0 v§ - 

p) = pe cos E 0 

(15c) 


The general form of Kepler 1 s equation and the general 
expressions for r, (rf) and (rv 2 -p) are obtained by 
substituting: 

r 0 , (r 0 r 0 >, (r 0 v 2 -y) Riven by (15a) - (l|c)» 


defined by (2) , 


11 


Sq , S| , S 2 , S, given by (6a) - (6d) and 
(2h) defined by (3), into equations (13) 
and (l^a) - (l*Jc). 


Kepler's equation for the general case is thus 


(t~t 0 ) * r 0 * + (r 0 r 0 ) S 2 + (r 0 vg - p)S 3 (16) 


in terms of \ P and its transcendental functions S 2 
and S 3 . Also, r, (rr), and (rv 2 - p) are 


r - r 0 + (r 0 r 0 )Sj + (r 0 vjj « p)S 2 ' (17a) 
(rr) * (r 0 r 0 )S 0 + (r 0 v§ - p)Sj (17b) 
(rv 2 - p) * (r 0 r 0 )(2h)S! + (rov* - w)S 0 (17c) 


in terms of- the transcendental functions S 0 , Sj , S 2 , 

S 3 . Use of the properties of (7a) shows that r, rr, 
and (rv 2 - p) are the successive derivatives of 
(t~t 0 ) with respect to ij> . This fact facilitates 
the graphical interpretation of Kenler's equation 
(16) in which (t-t 0 ) is plotted as a function of 
for all the various types of orbits. These Include 
the circular, elliptic, parabolic, hyperbolic, and 
rectilinear cases for an attractive force, as well 
as the hyperbolic and rectilinear cases for a repulsive 
force. 


12 


1.4, Closed Form Expressions for the Series Solution 

The series solution to the differential equations (1) 
is usually expressed in the form 


Cx, 

[x. 


y> 

z] =■ f*[ 

x o * y 0 » ?, o-^ 


V() 9 Z Q ] 

( 18a) 

y> 

z] - f-r 

*0> ?0» z 0^ + 

e-Cx 0 » 

o 

N • 
o 
UJ 

( 18 b) 


where f, g' and their time derivatives f, g are 

infinite power series in (t-to) whose successive 

coefficients are Increasingly complicated functions 

of x 0 , yo* z 0 , x 0 , y 0 > z 0 . However, the functions 
• • 

f, g, f, g can be expressed in closed form, in terms 
of the parameter (E-E 0 ) v for an elliptic orbit. 


f * 1 _ l-cos(E-E 0 ) 
r 0 7a 

fr = (t-t 0 ) - (E-Eq) - sin (E-Ep ) 

/^ 3 


-■J! 


g = l - 


sln(E-E n ) 
3 (r 0 /a) (r/a) 

1 - cos(E-E 0 ) 
r/a 


(19a) 

919b) 

(19c) 

(19d) 


These formulas are easily expressed in terms of S 1# 
S 2 , S 3 by use of the definitions (6b), (6c), (6d). 
The closed-form expressions for f, g, f, g in the 
general case are thus 
f = 1 - yS 2 /r 0 
g = (t-t 0 ) - pS 3 


( 20 ) 
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f - -ySj/(rr 0 ) 
g - 1 - wS 2 /r 

in terms of the transcendental functions Sj, S 2 , S 3 and 
the parameters r, r„ , y and (t-t 0 ). 
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2, A METHOD FOR COMPUTATION OF COORDINATES 
2.1. Initial Computations and Start of Iterations 

The first step In computing the numerical values of 
the coordinates Xj, y l9 z l9 Xj, y l9 i\ at a given 
time ti from the given coordinates x 0 , yp, z 0 , 
x 0 » Yo> z 0 at the given time t 0 is to calculate 
the parameters 


r 0 - + /x 0 =4 7 ,,*+ z“ 2 (21a) 

(r 0 r 0 ) ■ x 0 x 0 + y 0 y 0 ZqZq (21b) 

(vj) ■ x§ + yjj + z\ (21c) 

( 2h) » v* - 2y/r o ( 21d) 

(r 0 v „- H ) ■ r Q (vj) - u ( 2 le ) 

( t j — t q ) ■ tj — to (2If) 


from the coordinates x 0 , y 0 , z 0 , x 0 , v 0 , z 0 , the 
constant m , and the times t 0 and tj . 

When (2h) obtained in (21d) is negative and the orbit 
is thus periodic, the neriod 

P = 2/ry/(/^2h) 3 (22) 

of the orbit is calculated and 

m * INTEGER portion of [ (tj-t 0 )/t> + 1/2] (23) 

is determined. This minimizes the absolute value of 

(24) 


( t i - 1 0 ) * ( t j —t 0 ) -mo 
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which is calculated and used in the place of (tj-t 0 ) 
to determine Xj , y X9 z l9 y x$ z x in all the com- 
putations that follow. Because of this, no distinc- 
tion is made between ( t ^ — 1 0 ) and (tj-t fl ) or quan- 
tities computed from either in formulas given below 

• • c 

to determine x 1 , y x , z X9 x x , y x> z x . 

Kepler's equation (16) must be solved by an iterative 
scheme to determine the value » p x for ip which cor- 
responds to (t r -t 0 ). That is, that value ^ 
of fp must be found which makes the right hand side 
of Kepler’s equation (16) equal to ( t x — 1 0 ) • In 
describing these computations below, t/> will be used 
to represent a current approximation for i/ij , and 
\p 1 will be used to represent a further approximation 
calculated from the approximation ii> . The initial 
value for ^ is 

\p * (ti-t^J/pQ (25) 

which is computed and then used to evaluate the trans- 
cendental functions as described in 2.2. 

2.2. Evaluation of the Transcendental Functions 

The transcendental functions and S 5 are first 

computed from the current approximation ip for , 

by using equations (6e) and (6 f) in the forms 



i|i 1, [l/4! + (2ht|i 2 )/6! + (2hi|i 2 ) 2 /8l + (2hi|i 2 ) 3 /10! + 


S, - i)i 5 Cl/5! + (2hi)i 2 )/71 + (2hi|i 2 ) 2 /9! + ...] 


( 26 ) 


The accurate computation of each of the series in 
brackets in these two equations is an important nu- 
merical problem. A simole solution is to forward sum 
each of the series term by term until the addition of 
another term does not change either sum. Then the 
accuracy of the summations may be improved if desired 
by backward nesting the same number of terms used in 
the forward summations . Multiplication of the two sums 
by and t|> 5 respectively then gives S 4 and S 5 . 


The functions S 3 , S 2 , S 1? S 0 are then computed by 
using the relations 


S, = 


S 2 
Si 
S n 


i> 3 /6 + (2h)S 5 
<l> 2 / 2 + (2h)S 4 
it + (2h)S 3 
1 + ( 2 h ) S 2 


(27) 


which are obtained from equations ( 6f ) back through 
(6c). The functions S 2 and S 3 could be computed 
directly by two equations similar to (26) above, and 
could then be used to compute x t , y A , z 1 , Xj, y 19 z ^ . 
However, S4 and S 5 cannot in general be accurately 
computed from S 2 and S 3 , and S 4 and S s are 
required if partial derivatives are desired. The 
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functions S 0 and S x are defined and comouted 
merely for convenience of notation* 

2.3. The Solution of Kepler's Equation 

The value (t-t 0 ), corresponding to the value $ 
and its functions S 2 , S 3 , is first computed by 

(t-t 0 ) x r 0 tj> + (r 0 r 0 )S 2 t (r Q v2 - y)S 3 (28) 

which is Kepler’s equation (16). That is, if ( t — 1 0 ) 
were the time interval at which a solution were de- 
sired, \p would be the solution of Kepler’s equation. 
However, the iterative procedure must find that par- 
ticular ip for which the residual 

At * (t«t Q ) - ( t x — t Q ) (29) 

is zero. This particular value of ^ will then be 
the correct value for ^ . The residual (29) for 

the current value of if* is computed along with the 
current value of r which corresponds to ij» and 
its functions Sj, S 2 . 

! 

r ** r 0 + (r 0 r 0 )Sj + (r 0 VQ - y)S 2 (30) 

This r, is also the derivative d ( At ) /di|# and is 
therefore the slope of the curve of (t-tq) as a 
function of i|> . 
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Newton's method is now applied to determine a new 
approximation for ^ . 

i|>* - At/r (31) 

Then the transcendental functions SJ, S J and S J , 

SJ, S j, SJ of are computed by applying the 
formulas in. Section 2.2. but using i|>* rather than 
*P . Also, the results are used in equations (29), 

(30), and (31) above to obtain values (t * -t 0 ) , At', 
and r f which correspond to '!>' ' and its functions 
SJ, SJ S J . If the residual At* is then less in 
absolute magnitude than At, then ip f , the transcendental 
functions SJ, SJ, SJ. SJ, SJ, SJ and the functions 
(t ' -t 0 ) , At ' , r f are all accepted as new values for 
ip and S Q i S p S S 3 , S ^ , S g ( t *t q) } At, r . 

Then Newton’s method (31) is used to compute a nev; 
ip ' and repeat the entire computation. 

When At f is not less In absolute magnitude than 
At, the current ' is not accepted as a new value 
for \p . Rather, a different value for \p f is com- 
puted by setting n equal to unity in the equation 

ij>' • ip - At/n (32) 

Then the slope r has been replaced by unit slope 
to determine the new approximation ip* for ^ . 
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This Is then used to compute Its functions 

S(5, S», , S*, S» and (t»-t 0 >. At', r*. If At* 

is then, less In absolute magnitude than At, ^* 
and Its functions are accepted as new values for 
and Its functions, and Newton's method (31) Is again 

applied as described above* 

* 

* 

If At* Is not less In absolute magnitude than At, 
n In equation (32) is doubled to compute a different 
i|>* . The slope n is repeatedly doubled until a At* 
Is obtained which is less than At In absolute mag- 
nitude, or until and i|>* are numerically Identi- 
cal. When the latter Is true, the resulting Is 

accepted as the value of for computing the coor- 

. * 

dinates x lf y Jt z r , x x , z y . 

2.-4. Computation of Coordinates 

Since $ is now the correct value for ^ i , the func- 

• • 

tions f, g, f, g in equations (20) are the functions 
for the coordinates Xj , y 19 z x , Xj, y x , z x . There- 
fore, these functions are computed from (tj-tp) and 
the functions S 1# S 2 , S 3 and r x of ^ i 
* f « 1 - yS 2 /r 0 

g « (t x -t 0 ) - mS 3 
f * -pS 1 /(r 1 r 0 ) 

g « 1 - vSj/rj 


( 33 ) 
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*« *» 

2 # , ** 


These functions are then used to compute x,, yj, z,, 
xj, yi , Zj by using equations (18a) and (18b). 

x, - fx 0 + gx 0 

y, ■ fy 0 + S.v 0 

Zj - fz 9 + gz 0 (3 1 *) 

X! - fx 0 + gx 0 

v, - fy 0 + gy 0 

z, - fz 0 + gz 0 

Thus the coordinates x t » Vj > z,, x 1 , y t , Zj at 

time t | have been computed from the constant y , 

the time t 0 and the coordinates x 0 , v 0 , z 0 , x Q , 

• • 

y 0> z 0 at time to • 
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3. A METHOD FOR COMPUTING PARTIAL DERIVATIVES 

0 

3.1. Outline of Derivation of Partial Derivatives 

Let the times t 0 and t x as well as v be treated 

as fixed constants and let the coordinates x lt , z^, 
yj* zi be treated as dependent variables of x 0> 
y , z 0 , i fl , y Q , i 0 which are treated as independent 
variables. The thirty-six partial derivatives of each 
of the coordinates x p y t , z a , x , y 1 > z^ with respect 
to each of the coordinates x Q , y 0> z 0 , x Q , y Q , z Q 
have many important practical applications. These 
derivatives are obtained by chain differentiation of 
the relations in Section 2 that are used to compute 
y j » Zj, y j s f r om Xq 5 Y q » x q , y q » z o * 

The chain differentiations must then be combined to 
obtain tractable formulas for computing the partial 
derivatives . 

The chain differentiations are rather tedious and 
lengthy and will therefore not be given here. The 
whole procedure is facilited by the use of matrix 
notation. The basic idea is to obtain matrix rela- 
tionships between all the differentials of quantities 
which are direct or indirect functions of x 0 > Yo> z o> 
x 0> y 0 , z 0 « These results are then combined to 
eliminate all differentials other than dXj, dyj, dZj, 
dx,, dy j , dZj and dx 0 , dy 0 , dz 0 , dx 0 , dy 0 , dz 0 . Then 



- 22 


the coefficient matrix relating these differentials 
is the desired matrix of the thirty-six partial de- 
rivatives . 

3.2. Evaluation of Parameters and Periodicity Computations 
The parameter r x and the transcendental functions 
S 0 j Sj, S 2 , S 3 , S 4 , S g of ^ have been determined 
in the computations described in Section 2 to obtain 
the coordinates x 1 , , z x , x l , y ^ , i l . The para- 

meter rj must also be computed from 

r i = r( r oV s o + (r o v o -> 1 > s i3/ r i (35) 

as as 2 as 3 

and r-tr^j -r— , — must be computed from (7b) 
o n o n 3n 

1?JL * ^1^2 ”^3 f 

&h 

^ - * 1 S 3 - ( 36 ) 

= *jS 4 - 3S s 

In addition, the true values of ^ , S 3 4 S 4 , S g 
must be determined from ^ , S 3 , S 4 , S g by using the 
equations of (12) if the orbit is elliptic and m 
is not zero. In the computations of Section 2, 
no distinction was made between ^ , § 3 , . S 4 , S 5 
computed from C t x — t Q ) and the true values , S 3 , 

S 4 , S 5 which are the same functions of ( tj -t 0 ) . 
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However, this distinction must be made in order to 
obtain the correct partial derivatives from the 
formulas in Sections 3*3« and 

^ * iji *f m (2ir//-2h) 

S t * S 3 + m[2ir/(/-2h) 3 ] 

* S h + m[2ir/(/-2h) 3 ](if» + I?)/2 (37) 

S 5 * S 5 + m[2ir/(/-2h) 3 ][(!? 2 +^+^ 2 )/6 + l/2h] 


The functions S 0 , S 1$ S 2 need not be recomputed 
since they are equal to S 0 , Si, S 2 by (12). 

3.3. Evaluation of the T?our by Four Matrix 

The four by four matrix below is first calculated as 
an intermediate step for the computation of the par- 
tial derivatives. The letter T is used to indicate 
the matrix transpose of the two column matrices. 


a l 1 

a i2 

a l 3 

a l 4 

a 2 1 

a 2 2 

a 2 3 

a 24 

a 3 1 

a 3 2 

a 3 3 

a 34 

a 4 1 

a 42 

a 4 3 

a 44 


0 

0 

s l /r l. 

V r 2 


( S q - r 1 S j ) /r 1 
(S j -fjS 2 )/rj 


h +V 0 S 3 + ?2 f( ^0 ^0 > ( r 0 V 0 -V 

ft 


r 0 S 2 


_ _ r . , 3S,.. 2 .9S,1 

2i , o s 3 + [(r 0 i'o>3h L+ (r 0 v 0 -w) 517 iJ 


. T 


1 +V 0 S 2 + ^ + (r 0 V 0-K)’3-^- 


r 0 S l 


r 0 S 2 + 


[(r 0 r 0 )||i + (r 0 Vg - v ) 


9S 

9h 


] 
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W / V 3jS 

r 0 \r 0 ah 



0 



0 


0 

0 




0 


V * S i 

*T i TK" 


(n/ r o) 2 as 2 

— 5TT 



y jS 

FT ih 


3.4. Computation of the Partial Derivatives 

The two by two sub-matrices of the four by four matrix 


above are then used to calculate the partial deriva- 


tives by the formulas given below. The letter T 
indicates the matrix transpose of the three by two 
matrices . 


3x 3x 3x 
*x 0 3 y 0 dz o 



*CL 

ro 


3y 3y 3y 

HliMlWI ««M|M 

3x„ 3y 0 3z„ 



3 Z 3 2 3 Z 

3x 0 3y 0 3z 0 





z 

— 0. z 
ro o 




< 

■i 

j 

4 

% 

j 

i 


1 

A 

'l 

f 



3x 3x 3x 

7x7 »y 0 9 


5y_ ii_ is_ 

9X„ 9y 0 9*0 

3z 3 Z 3 

9*0 9y 0 9Z 0 . 


3x 9 X 3x 

^ *77 ^ 


**7 

*77 

*77 

3z 

3z 

3 Z 

^ 57 ; ^7 

3X 

3 X 

3x 

x 0 

3 yo 

3Z 0 

ii.. 

ay 

ii_ 

• 

3X 0 

3yo 

3Z 0 

3z 

• 

3 z 
• 

3z 

• 


3x« 3v* 3 z 
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CONCLUDING REMARKS 

The authors have written a FORTRAN IV program (available 
upon request) for computing coordinates and partial deri- 
vatives of the two-body problem. Cases run to test the 
program included: (for v>0) elliptic circular, parabolic, 

hyperbolic, rectilinear; and (for wiO) hyperbolic, recti- 
linear. 

Computationally, the program is superior to available 
programs in that it produces solutions and partial de- 
rivatives for all cases of the two-body problem without 
exception. It also has no disadvantage in the accuracy 
and speed of computation. 
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